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Abstract 

We study quantum oscillations of the magnetization in Bi2Se3(lll) surface system in the presence 
of a perpendicular magnetic field. The combined spin-chiral Dirac cone and Landau quantization 
produce profound effects on the magnetization properties that are fundamentally different from 
those in the conventional semiconductor two-dimensional electron gas. In particular, we show that 
the oscillating center in the magnetization chooses to pick up positive or negative values depending 
on whether the zero-mode Landau level is occupied or empty. An intuitive analysis of these new 
features is given and the subsequent effects on the magnetic susceptibility and Hall conductance 
are also discussed. 
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Magnetic oscillation, which was first predicted by Landau in 1930 jlj|, has been a focus in 
the condensed matter physics filed. One important reason is that the de Haas-van Alphen 
(dHvA) oscillations of magnetization provide a vigorous technique to study the properties 
of carriers around the Fermi surface. Especially in the last decade, thanks to the tremen- 
dous advances in microscopic semiconductor technology, the challenge encountered in the 
measurement of weak magnetization signal has been largely prominently overcome, and the 
magnetic dHvA oscillations in the two-dimensional electron gas (2DEG) systems have thus 
been extensively studied. For instance, Meinel et al. [2H^ developed dc superconducting 
quantum interference device magnetometers to study the dHvA oscillations in high-mobility 
semiconductor 2DEG. Schwarz et al. studied the dHvA oscillations by using mi- 

cro mechanical cantilever magnetometers. Besides the purely orbital part, prominently, the 



influence 



of the weak Rashba spin-orbit interaction (SOI) on the dHvA oscillations in 



the magnetization of the semiconductor 2DEG can also be effectively determined in experi- 
ment [91 , which therefore opens a new door to measurement of the spintronic parameters in 
semiconductor heterostructures. 

In the above-mentioned conventional semiconductor 2DEG systems, in which the elec- 
tron motion is dominated by its orbital part, i.e., the magnetization oscillation mainly comes 
from the response of the electron /c-quadratic kinetic energy to the external magnetic field. 
Although sometimes other factors than the pure kinetic energy, such as the spin-orbit inter- 
action Sj, have been taken into account, these factors in conventional semiconductor 2DEG 
systems play only a minor role. For example, they can result in a beating mode superposed 



onto the main dHvA oscillation pattern \0^. This situation, however, wi 



1 be greatly changed 



111 ] of the topologi- 



by very recent theoretical finding [lO| and experimental verification 
cal insulators (TIs) with strong spin-orbit interaction. As a new state of matter as first 
addressed by Kane and Mele [l2|, the subject of time- reversal invariant TIs has attracted 
great attention in condensed-matter physics. Several three-dimensional (3D) solids, such as 
Bii_a;Sba; alloys, Bi2Se3-family crystals, have been identified |l3l-ll7| to be strong TIs possess- 
ing anomalous band structures characterized by a Z2-valued topological invariant 12, Il8 |. 
This invariant, called z/q, counts the number of topologically protected gapless surface states 
(modulo 2). A non-zero invariant means that the surface of 3D TIs will be metallic. Instead 
of the conventional semiconductor 2DEG systems, the energy scale for the surface states 
of these 3D TIs is dominated by the fc-linear spin-orbit interaction instead of the parabolic 
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FIG. 1: (Color online) (a) The ab initio calculated band structure of the Bi2Se3(lll) thin film 
with the thickness of six quintuple layers. The red lines indicate the surface states while the black 
lines correspond to the bulk bands, (b) The fitting curves (blue lines) around the T point with the 
model Hamiltonian [Eq. (1)]. 

kinetic energy. As a result, it is expected that the magnetic response properties of these 
topological surface states are fundamentally different from those of the conventional 2DEG. 
Inspired by this observation, as well as by the recent experimental observation of the 

19|, in this paper we study the elec- 



Landau quantization of the surface states of Bi2Se3 
tron magnetic oscillations of these surface states. Specially, we consider the surface states 
of Bi2Se3. The first-principles surface band structure of Bi2Se3 is calculated by a sim- 
ple supercell approach with spin-orbit coupling included and shown in Fig. [U^a) along the 
high-symmetry lines (F— t-M, M— j-K, and K— J-F) in the surface Brillouin zone. In obtain- 
ing Fig. [U^a), here we have used Vienna ab-initio simulation package (VASP) [20 1. The 



(Perdew-Burke-Ernzerhof) 
augmented wave potential 



'BE [2X] generalized gradient approximation and the projector- 
22[ were employed to describe the exchange-correlation energy 
and the electron-ion interaction, respectively. The SOI, which has been confirmed to play 
an important role in the electronic structure of Bi2Se3, was included during the calculation. 



3 



The cutoff energy for the plane wave expansion was set to 300 eV. During the calculation, 



the experimental lattice constants 



23| were adopted, i.e., a=4. 143 A, c=28.636 A, with 



the internal parameter optimized automatically. The Bi2Se3(lll) surface was modeled by 
a slab composing of six quintuple layers (QLs) and a vacuum region of 20 A. Integration 
over the Brillouin zone was done using the Monkhorst-Pack scheme 2J] with 10x10x1 grid 
points for surface calculations. The structures of bulk and slab were fully optimized until 
the maximum residual ionic force were below 0.01 and 0.02 eV/A, respectively. From Fig. 
[U^a) two chiral surface states are clearly seen to connect the conduction band and valence 
band, and cross each other to form a single Dirac-type contact at the F point and aligning 
with the Fermi energy. The intrinsic defects such as the substitutional Bi defects at Se sites 
or the Se vacancies will play a role of ra-doping, and consequently shift the Fermi level above 
the Dirac point. Note that due to the difference in the local symmetry between the top and 
bottom surfaces, there can develop an observable asymmetric charge distribution on the two 
surfaces if the sample is thin enough. This fact sometimes can open a small gap between 
the two spin chiral surface states as shown in Fig. [I](b) which presents an enlarged version 
of the surface bands around the Dirac cone. With increase of the thickness of the film, 
however, this asymmetry-induced Dirac gap tends to vanish, and this actually corresponds 
to the recent Landau quantization experiment, in which the used epitaxial Bi2Se3 is as thick 
as 50 QLs. 

We report the calculated magnetization of the electrons on Bi2Se3 surface as a response 
to the external magnetic field. It is found that the magnetization oscillations in the present 
system differs from the traditional 2DEG by the fact that the dHvA oscillating center in 
Bi2Se3 departs from the well known (zero) value in the semiconductor 2DEG system. This 
departure has an intimate relation with the different Landau spectrum structures in these 
two kinds of systems. It is well known that the Landau levels (LL's) in the traditional 
2DEG obey the B{n+l/2) rule with B being the external magnetic field and n the LL 
index. However, the energy spectrum of the surface states in Bi2Se3 approximately obeys a 
y/nB rule. It is this difference in LLs that distributes differently in the two components of 
the magnetization, and eventually result in the different magnetic properties in these two 
systems. Furthermore, we have shown that the zero-mode LL plays a key role in determining 
the magnetization behavior in the TI surface systems. 

The Hamiltonian describing the gapless surface states of Bi2Se3 can be approximately 
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written as 

H{k) = + hvp [kxCTy - kyax) , (1) 
where is the Fermi velocity and a are the Pauh matrices for surface-state electron spins. 



Note that this Hamiltonian has t 
with Rashba spin-orbit coupling 



le same form as that of the conventional 2DEG system 
8|. However, the intrinsic difference between these two 
kinds of systems is that the fc-linear spin-orbit interaction is primary to the TI surface 
states, while the parabolic term is dominant in the conventional 2DEG. Although it is very 
simple, the Hamiltonian ([2]) is a general one, which can satisfactorily describe the gapless 
surface states of Bi2Se3 near the Dirac point. This satisfaction is particularly obvious for the 
upper part of the Dirac cone (electron part), as shown in Fig. 1(b), where Eq. ([2]) is used 
to fit the first-principles result, which gives 7=2.1x10^ meV nm^ and ^fir=200 meV nm 
(namely, ^^=3.04 x 10^ m/s). The effective mass m* is then obtained as O.lSme according 
to 7=/i^/2m*, where rrie is the mass of a free electron. The lower surface states (hole part) 
is not well described by Eq. ([2]) and a better fitting needs higher fc-order corrections. For 
simplicity, and for the reason that we only concern the n-doping, here we neglect 0{k^) 
terms. 

Let us now consider an external magnetic field B=Bz being perpendicular to the surface. 
Taking the Landau gauge for the vector potential, Ax=By and Ay=0, and the transform 
fik-^Il=fik+eA, one can obtain the following Hamiltonian 

n2 1 

H = ^— ^ + vf {Jlx(yy - Llj^cr^) - -gs/J'sBa^, (2) 
where Qs is the effective magnetic factor of the surface electron and /i^ is the Bohr magneton. 



For Bi2Se3-family (111) surfaces, the value of gs is approximately 8.0 [25|]. Taking the fact 
that the system is translation invariant along the x axis and therefore the wave number 
along this direction is a good quantum number, the Hamiltonian ([2]) can be rewritten as 

H = tkUc ^a^a H + i\Plr]{aa^ — aV+)^ , 

where (j±=((Ta; ±icry)/2, ujc=eB /m* , ri=VFm*lB/fi is the effective Rashba spin-orbit coupling 
with lB=\Jfi/eB being the magnetic length, g=gsm* /2me, and a=[y+{hkx + ipy)/eB]/y/2lB 
is the usual harmonic oscillator operator. The LLs are then given by 

= (^n ± i V(l-(7)2 + 8V^ (3) 



with 72=1, 2, ■ ■ ■ . The n=0 LL only has the "+" branch, Eq'^^ =^-^-^HuJc- The corresponding 
two-component eigenstates for E^'' are given by 

( / cos^i^V) \ 
\ 2 sin 6'n p — 1) / 

where \n) is the eigenstate of the nth LL of a free two-dimensional spinless electron. Here, 
tan ^1"^''=— u„ ± + with — g)/V8nri when n>0 and ^q'''''=0 for ?2=0. It is 

interesting to see that the n=0 LL has the fully polarized spin along the z direction. Figure 
[2] plots the LLs as functions of the magnetic field. Note that although the LL equation 
([3]) for Bi2Se3 surface states has a similar form with the conventional spin-orbit coupled 
semiconductor 2DEG sl, these two systems are fundamentally different by the amplitudes 
of the physical parameters. For the former the dimensionless parameter rj^l while for the 
latter r^^l. Actually, for Bi2Se3, the Se-terminated (111) surface lattice constant is a=4.143 
A. With this knowledge and through a normal fitting process, we obtain that at the external 
magnetic field B=l T, fluJc=0.61 meV and ?7=12.4. However, for a conventional 2DEG 
system with Rashba coupling, the dimensionless parameter rj is typically in the range 0~0.2. 
Based on this fact, the energy spectrum ([3]) can be well approximated by the dispersion 
relation 



E^^) = ±noJc^ 2nr]^ + C " ±VFV2nehB + 6^ {n ^ 0), (5) 
4+) = sgn(^?,K|5|, 

where 5=gs^BB /2vf- Because the Zeeman splitting is much smaller than the LL separations 
(for example, (^=0.72 when Qs takes the value 8, resulting in |(yfs/iBi?=0.13 meV at B=l T), 
thus the effect of the Zeeman term on the n^Q LLs is very tiny and can be safely neglected 
in considering the electron occupation of n^^O LLs. It is not so, however, for the ?t,=0 LL. 
In fact, in the absence of the Zeeman splitting, the Dirac-Landau energy spectrum ([5]) is 
massless with a whole electron-hole symmetry and only half of the zero mode is occupied 
by electrons in the case of n-doping. If the Zeeman splitting is finite, the spectrum ^ 
is massive and the n=Q LL shifts upward or downward, depending on the orientation of 
the exchange field (the sign of Qs)- Correspondingly, this "zero" mode will be saturated 
by electrons for gs>^ or empty for 5's<0, which, as shown in the following discussion, will 
greatly influence the behavior of the magnetic response of the system. 
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FIG. 2: The Landau levels in the Bi2Se3(lll) surface system as functions of the external magnetic 



Now with the LL spectrum ([5]) [or the fc^-corrected LL spectrum ([3])], we study the 
magnetization of the surface electrons of Bi2Se3. The magnetization density is the derivative 
of the Helmholtz free energy density with respect to B at fixed electron density A/" and 
temperature T, M=—{dF/dB)\j^^T- The Helmholtz free energy density is given by 



where [5=1/ UbT, Ny=l/2Til\ is the degeneracy for a non-zero LL (namely, the number of 
electrons per unit area on a LL), and /i is the electron chemical potential. The second line 
in Eq. ([6]) denotes the contribution from the n=0 LL and there exist three possibilities for 
its contribution: (i) If this level is exactly a zero mode, i?g^''=0, then the system has the 
electron-hole symmetry and half of the particles in the zero mode are electrons. In this case 
Nq=N^/2; (ii) If the Zeeman splitting cannot be neglected and the g factor is positive as Eq. 
(IH]) shows, then the n=0 LL shifts upward and is fully accessible to electrons. In this case 



field B. 




n=l 





Nq=N^; (iii) Otherwise, if the g factor is negative, then the n=0 LL shifts downward and is 
unavailable to electron occupation. In this case A^o=0. The -B-dependent chemical potential 
fi in Eq. (E]) is connected to the experimentally accessible electron density M, which is given 
by 

oo 



n=l 



with /^+^=1/ 



1 



being the Fermi-Dirac distribution of the LL Er^\ From Eq. 



the electron magnetization density becomes 




OB 




(3 dB 



1 + 



(8) 



The first part M^^^ is the conventional contribution from the B dependence of the LLs and 
thus denotes a diamagnetic response. The second part M^^^ comes from the B dependence 
of the level degeneracy factor N^, thus describing the effect of the variation in the density 
of states upon the magnetic field and denoting a paramagnetic contribution to the total 
magnetization. Obviously, M*^°) is negative while M*^^) is positive, the net result is an 
oscillation of the total magnetization M as a function of B. 

We plot in Figs. |3]^a) and |3]J^b) the magnetic dHvA oscillations of the chemical poten- 
tial and magnetization in Bi2Se3 for the above-mentioned three cases of zero-mode filling. 
Comparing to the well-known dHvA oscillating pattern in the conventional 2DEG, one can 
immediately obtain three prominent features in the present TI surface system: (i) Although 
the oscillating center of the chemical potential /x keeps a fixed value unchanged by chang- 
ing the external magnetic field B when the n={] LL is exactly a zero mode, it linearly 
increases/decreases with B when the g factor is positve/negative; (ii) The oscillating center 
of the magnetization M keeps a non-zero value unchanged by varying the external magnetic 
field strength. This is totally different from those in the semiconductor 2DEG. It is well 
known that in the clean semiconductor 2DEG sample, the oscillating center of the chemi- 
cal potential ^ keeps a constant value unchanged and that of the magnetization M keeps 
zero unchanged when varying the magnetic field B] (iii) The magnetization for the case of 
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FIG. 3: (Color online) Magnetic field dependence of (a) chemical potential and (b) magneti- 
zation M in n-doped Bi2Se3(lll) surface system with electron density A/'=1.0 x 10^^ m~^. The 
temperature is set as kBT=0.3 meV. The black, red, and blue curves correspond to the cases that 
the zero mode is half filled, saturated, and empty, respectively. 



empty zero mode is fundamentally distinguished from the cases of saturated and half-filling 
zero mode by a total sign inversion. In addition, the magnetic oscillation patterns for the 
cases of saturated and half-filling zero mode are out phase. Thus, Fig. [3] clearly reveals the 
fundamental role the zero mode played in determining the magnetic response properties of 
the TI surface states. 

For further illustration and to see the origin of the sign inversion in the magnetization 
when the zero is changed from filling to unfilling, let us first consider the case of saturated 
zero mode. In this case the n=0 LL is occupied by electrons with the degeneracy Nq=N^. 
The discrepancy in the dHvA oscillating patterns between the TI surface and the semi- 
conductor 2DEG comes from their different energy dispersion relations. The former versus 
the external magnetic field obeys square root rule while the latter obeys linear rule. It is 
well known that at zero temperature, the two components of the magnetization in the con- 
ventional 2DEG turn to be M(o)=-f J^lZo K and M(i)=f En=o'(/^o - E„), respectively. 
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Here /^o is the zero-temperature Fermi energy and the LL En ccB. The negative M^°^ and 
the positive M'^^^ gives that the net result is an oscillation of the total magnetization M 
as a function of B. The oscillation amplitude increases with B and the oscillation center 
is zero. However, because the LL EnOcy/B for the present system, at zero temperature 
the first component of M turns to be M ^ X]n=o ^2~' *hile the second component is 
^^^■' = f Sn=o (/^o ~ E^^). By comparison with those in the semiconductor 2DEG, one can 
find that in the present TI surface system the diamagnetic contribution (M'^°)) is reduced. 
As a result, the oscillating center of the magnetization is now a positive value for the case of 
saturated zero mode. This simple comparison is not strict in mathematics, however, it af- 
fords an intuitive explanation on the difference of the magnetization between the TI surface 
and the semiconductor 2DEG. 

In the case of empty zero mode, the n=0 LL is excluded and the first component M*^°) now 
becomes M'^°^=— | while the second component becomes M*^^) = | ^°^"'(/io — 

E^'^). Compared to the case of saturated zero mode, and considering ijlq':^e'^\ one can 
find that the magnetization in the case of the empty zero mode is smaller than the saturated 
case by a quantity |/io and therefore becomes negative during its oscillations as a function 
of S. 

Note that the abrupt jump in the dHvA oscillation is on the high magnetic field side 
of the sawtooth, which is special for our present choice of the thermodynamic system. If 
the system is constrained to have constant chemical potential, then the jump in the dHvA 
oscillation will move to the low magnetic field side of the sawtooth, which has been confirmed 
by Meinel et al. in an experiment with the electron density M modulated by applying a 
gate voltage to the sample. 

The above discussions on the dHvA oscillations of the magnetization focus on the situa- 
tion that the total number of electrons on the LL's is field independent. Now we consider 
the magnetization properties in another situation that the chemical potential is field inde- 
pendent. Figure 111(a) plots the magnetization as a function of the zero-temperature chemical 
potential (i.e., the Fermi energy) at 5=4 T for three cases of zero- mode filling. The dHvA 
oscillating patterns as a function of the chemical potential can be observed from Fig. Il](a). 
There exist different kinds of patterns for the dHvA oscillating center in different cases. 
When the n=0 LL is half occupied by electrons, the dHvA oscillating center keeps zero 
unchanged by increasing the chemical potential. When the zero mode is saturated/empty, 
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FIG. 4: (Color online) (a) Calculated magnetization as a function of the Fermi energy ^uq- The 
external magnetic field is chosen as B=4 T. (b) The derived Hall conductance from dM/dfiQ. The 
black, red, and blue curves correspond to the cases that the zero mode is half filled, saturated, and 
empty, respectively. 

the dHvA oscillating center linearly increases/decreases by increasing the chemical poten- 
tial. A fact that should be addressed is that because the chemical potential is tuned freely 
and independent with the external field, there are no phase difference in the dHvA oscil- 
lations for different cases. The corresponding dM/dfiQ are also calculated, from which the 
Hall conductance an is obtained by combining the thermodynamic relationship and Streda 
formula: 

^M/^^Io = -an. (9) 

e 

The result of Hall conductance as a function of the Fermi energy is plotted in Fig. IH^b), 
from which Hall plateaus can be clearly seen. The plateau values of an depend on the 
zero-mode filling. If the n=0 LL is half filled, the Hall conductance takes half-integer values 
of (jH={n+l/2)e^ /h, as shown in Fig. 111(b) by black step lines. To date, measuring the 
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half-integer quantum Hall effect on the TI surfaces keeps a challenging task, although the 



LLs have been recently observed 19|. If the zero mode is saturated, then aH={n+l)e^ /h, as 
shown in Fig. lU^b) by red step lines. Finally, if the zero mode is empty, then a H={n—l)e^ / h, 
as shown in Fig. IH^b) by blue step lines. We note that the quantum Hall effect in the TI 
surface system with finite sample size has also been discussed in Ref. 26 1. 

The information on the magnetic susceptibility x of the TI surfaces, which is defined as the 
derivative of the magnetization with respect to the external magnetic field, x=dM/dB, can 
be easily obtained with the knowledge of the magnetization [Eq. (|8])]. The final expression 
of the magnetic susceptibility is given by 



X 



/■(+) 




dB dB 



(+)" 



(10) 



dB dB 



Jo 



dB^^ 



Figure [5] plots the magnetic susceptibility in Bi2Se3 sample as a function of the inverse 
magnetic field, 1/B. Because the resonant (for magnetic susceptibility) magnetic field only 
reflects the occupation of the LLs, which is the same as that in the conventional 2DEG, 
it losses the message on the oscillating center value of the magnetization. However, the 
difference between the half-filled case and saturated/empty case of the zero mode is still 
clearly revealed in this figure. 

In summary, we have investigated the dHvA oscillations of the magnetization in the 
Bi2Se3-family surface systems. Our results show that the dHvA oscillating center of the 
magnetization maintains positive values when the n=0 LL is fully occupied or half occu- 
pied. When this mode is empty, the dHvA oscillating center changes a sign. These results 
are fundamentally different from those in the conventional semiconductor 2DEG systems, 
in which the dHvA oscillating center is at zero. We have given an intuitive analysis on 
this difference, which turns to have an intimate relation with different forms of the energy 
dispersions in these two kinds of systems. This can be reflected by the fact that the LLs 
for the TI surfaces is proportional to Vb instead of 5-linear dependence accommodated 
by the conventional 2DEG. Furthermore, the essential role that the zero mode played has 
been illustrated by the different behavior of the Hall conductance at three kinds of electron 
occupation of this mode. We expect that the present results for the topologically nontrivial 
features in the magnetic response of the TI surfaces can be experimentally confirmed in the 
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FIG. 5: (Color online) Calculated magnetic susceptibility x as a function of 1/B. The parameters 
are the same as those in Fig. [3j The black, red, and blue curves correspond to the cases that the 
zero mode is half filled, saturated, and empty, respectively. 

future topological magnetoelectric studies. 

Note added. — While this work was completed, we were aware of an experimental mea- 
surement [^^l of the magnetization for the topological insulator Bii„^.Sb^ (0.07<x<0.22). 
Compared to Bi2Se3, the surface band structure of Bii„a.Sba; alloy is much more compli- 
cated. Furthermore, in Bii-ajSb^ the bulk band is often coupled with surface band during 
measurement. These facts make it difficult to study the magnetic quantum oscillations that 
are fully from the 2D surface states of Bii-ajSb^,. alloy. In spite of these complicated facts, we 
expect that the phenomenon of large-amplitude dHvA magnetic oscillations found in Ref. 



27| is closely related to our theoretical finding in the present paper. 
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